In [1]:
# Get matplotlib magic
%matplotlib notebook
# Import modules
import numpy as np
from scipy import signal, interpolate
import matplotlib.pyplot as plt
from math import ceil
plt.rc('font', size=8)
# Run data extraction and useful function script; double quotes work for Windows and Mac paths
%run "data_and_funcs.ipynb"
# AJG functions
# help(get_erp)
# help(low_pass_erp)
# help(band_pass_erp)
Figure 9.1¶
In [2]:
# Specify channel label
which_chan = 'FCz'
# Find the index (channel number) of that label
chan_index = [index for index, item in enumerate(chan_names) if item == which_chan]
# Set up the figure and grid
fig_91 = plt.figure(figsize=(8, 4))
heights = [2.25, 1]
widths = [1.5, 1]
full_grid = fig_91.add_gridspec(nrows=2, ncols=2, height_ratios=heights, width_ratios=widths)
# Figure 9.1A ----------
# Determine the number of trials to plot, the total number of trials, and randomly select from available trials
num_trials_to_plot = 12
num_trials = eeg_data.shape[2]
random_trials_to_plot = np.random.permutation(num_trials)[0:12]
# Uncomment the following line if you would like to replicate Cohen's trials
# random_trials_to_plot = np.array([68, 74, 65, 70, 28, 10, 69, 95, 44, 76, 19, 45]) - 1
# Compute appropriate grid size for subplots
n_rows = ceil(num_trials_to_plot/ceil(np.sqrt(num_trials_to_plot)))
n_cols = ceil(np.sqrt(num_trials_to_plot))
# Set subgrid
grid_A = full_grid[0:, 0].subgridspec(n_rows, n_cols)
# Loop through random trials and plot them on the subgrid
j = 0
for i in random_trials_to_plot:
ax_A = fig_91.add_subplot(grid_A[j])
ax_A.plot(eeg_time, eeg_data[chan_index, :, i][0], color='black')
ax_A.set_xlim([-200, 1000])
ax_A.set_xticks([0, 500])
ax_A.set_yticks([])
ax_A.set_title('Trial ' + str(i + 1), loc='left')
plt.tight_layout()
# Label the bottom left corner subplot
if j == (n_rows * n_cols - n_cols):
ax_A.set_xlabel('Time (ms)')
ax_A.set_ylabel(r'Voltage ($\mu$V)')
j += 1
# Figure 9.1B ----------
ax_B = fig_91.add_subplot(full_grid[0, 1])
ax_B.plot(eeg_time, np.squeeze(eeg_data[chan_index, :, :]), color='gray')
ax_B.plot(eeg_time, get_erp(eeg_data, which_chan), color='black', linewidth=2)
ax_B.set_xlim([-300, 1000])
ax_B.set_ylim([-60, 60])
ax_B.invert_yaxis()
ax_B.set_title('All trials and trial-average (ERP)')
# Figure 9.1C ----------
ax_C = fig_91.add_subplot(full_grid[1, 1])
ax_C.plot(eeg_time, get_erp(eeg_data, which_chan), color='black')
ax_C.axhline(y=0, color='grey')
ax_C.axvline(x=0, color='grey', linestyle='dashed')
ax_C.set_xlim([-300, 1000])
ax_C.invert_yaxis()
ax_C.set_yticks([6, 4, 2, 0, -2])
ax_C.set_ylim([7, -3])
ax_C.set_xlabel('Time (ms)')
ax_C.set_ylabel(r'Voltage ($\mu$V)')
ax_C.set_title('ERP', loc='left')
fig_91.tight_layout(h_pad=2.5)
Figure 9.2¶
In [3]:
# Specify channel label
which_chan = 'P7'
erp = get_erp(eeg_data, which_chan)
nyquist = info['sfreq'] / 2
transition_width = 0.15 # percentage
# Run a low-pass filter from 0 to 40 Hz
erp_0_to_40 = low_pass_erp(erp, filter_cutoff=40, trans_width=transition_width, nyq=nyquist)
# Run a low-pass filter 0 to 10 Hz
erp_0_to_10 = low_pass_erp(erp, filter_cutoff=10, trans_width=transition_width, nyq=nyquist)
# Run a band-pass filter from 5 to 15 Hz
erp_5_to_15 = band_pass_erp(erp, filter_low=5, filter_high=15, trans_width=transition_width, nyq=nyquist)
# Figure 9.2 ----------
fig, ax = plt.subplots()
ax.plot(eeg_time, erp, linestyle='solid', color='black')
ax.plot(eeg_time, erp_0_to_40, linewidth=2, color='blue')
ax.plot(eeg_time, erp_0_to_10, linestyle=':', color='red')
ax.plot(eeg_time, erp_5_to_15, linestyle='--', color='green')
ax.set_xlim([-200, 1200])
ax.set_xticks([0, 200, 400, 600, 800, 1000])
ax.invert_yaxis()
ax.set_yticks([8, 6, 4, 2, 0, -2, -4])
ax.set_xlabel('Time (ms)')
ax.set_ylabel(r'Voltage ($\mu$V)')
ax.set_title('ERP from electrode ' + which_chan)
ax.legend(loc='upper right', labels=('None', '0-40', '0-10', '5-15'), frameon=False)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.show()
Figure 9.3¶
In [4]:
# Prepare the subplots
fig, (ax1, ax2) = plt.subplots(nrows=2, sharex=True)
# Axis 1
for i in np.arange(eeg_data.shape[0]):
ax1.plot(eeg_time, np.mean(eeg_data[i, :, :], 1), color="gray")
ax1.invert_yaxis()
ax1.set_ylim([15, -5])
ax1.set_ylabel(r'$\mu$V')
ax1.set_title('ERP from all electrodes', loc='left')
# Axis 2
ax2.plot(eeg_time, np.var(np.mean(eeg_data, 2), 0), color="gray")
ax2.set_xlim([-200, 1000])
ax2.set_yticks([0, 5, 10, 15])
ax2.set_xlabel('Time (ms)')
ax2.set_ylabel(r'var($\mu$V)')
ax2.set_title('Topographical variance from all electrodes', loc='left')
plt.show()
Figure 9.4¶
In [5]:
# We would just like to plot the topographical map from Figure 9.4
# Specify parameters
time_point_to_plot = 100
color_limit = 20
abs_times = np.abs(eeg_time - time_point_to_plot)
min_idx = np.where(abs_times == abs_times.min())
# Random trial; This is why it will not look exactly like the text book
trial_to_plot = np.random.randint(eeg_data.shape[2] + 1)
topo_data = np.squeeze(eeg_data[:, min_idx, trial_to_plot])
head_dict = {'center': [0, 0],
'scale': [.4, .4]}
# Use MNE to generate a scalp topography
plt.figure()
mne.viz.plot_topomap(topo_data, info,
vmin=-color_limit, vmax=color_limit, cmap='jet',
outlines='skirt',
contours = 8,
head_pos = head_dict
)
plt.show()
Figure 9.5¶
In [6]:
# Find the times in the EEG closest to the numbers in the sequence from -100 to 600 with step size of 50
times_to_plot = []
for i in np.arange(-100, 650, 50):
times_to_plot.append(np.abs(eeg_time - i).argmin())
# Set up the grid
fig_95 = plt.figure()
full_grid = fig_95.add_gridspec(nrows=3, ncols=5)
# Find the index for FC4
fc4_index = [index for index, item in enumerate(chan_names) if item == 'FC4']
# Iterate through the different times and plot the data averaged over trials
for i in np.arange(len(times_to_plot)):
eeg_data_plot = np.squeeze(np.mean(eeg_data[:, times_to_plot[i], :], 1))
eeg_data_plot[fc4_index] = np.random.normal() * 10
ax = fig_95.add_subplot(full_grid[i])
mne.viz.plot_topomap(eeg_data_plot, info,
vmin=-8, vmax=8, cmap='jet',
outlines='skirt',
contours = 8,
head_pos = head_dict
)
ax.set_title(np.array2string(np.arange(-100, 650, 50)[i]) + ' ms')
Figure 9.6¶
In [7]:
# Initialize reaction time vector
rts = np.zeros(eeg_data.shape[2])
for i in np.arange(len(rts)):
# Get the reaction times
time_0_ev = np.where(np.concatenate(epochs[i]['eventlatency'] == 0))
rts[i] = np.concatenate(epochs[i]['eventlatency'])[time_0_ev[0] + 1]
fig, ax = plt.subplots(ncols=2, figsize=[8, 4])
for i in np.arange(2):
if i == 0:
# Sort by reaction time
sort_idx = np.argsort(rts)
ax[i].plot(rts[sort_idx], np.arange(len(rts)), color='black', linewidth=2)
ax[i].set_title('ERP image sorted by reaction time')
else:
# Sort by voltage at 300 ms
sort_idx = np.argsort(np.squeeze(eeg_data[chan_index, 300, :]))
ax[i].set_title('ERP image sorted by voltage at 300 ms')
image = ax[i].imshow(np.squeeze(eeg_data[chan_index, :, sort_idx]),
extent=[np.min(eeg_time), np.max(eeg_time), len(rts), 0],
aspect='auto', vmin=-30, vmax=30)
ax[i].set_xlim([-200, 1200])
ax[i].set_xticks(np.linspace(0, 1000, 6))
ax[i].invert_yaxis()
ax[i].set_xlabel('Time (ms)')
ax[i].set_ylabel('Trials')
fig.colorbar(image, ax=ax[0])
plt.show()
